# -*- coding: utf-8 -*-
"""
Created on Tue Aug 31 13:05:22 2021

@author: zhiyinan
"""

import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import Parameters_1D as par_RNN 
from Polar1D_vectorized import *
import van_Genuchten1D_vectorized as vg 
from scipy import integrate
import matplotlib.font_manager as font_manager
from sklearn.preprocessing import MinMaxScaler
import random
from casadi import *
import tensorflow.keras.backend as kb


# %% Define NN layers
def sigmoid(x):
    return 1/(1+np.exp(-x))

def LSTMlayer_sig(warr, uarr, barr, x_t,h_tm1,c_tm1):
    
    '''
    c_tm1 = np.array([0,0]).reshape(1,2)
    h_tm1 = np.array([0,0]).reshape(1,2)
    x_t   = np.array([1]).reshape(1,1)
    
    warr.shape = (nfeature,hunits*4)
    uarr.shape = (hunits,hunits*4)
    barr.shape = (hunits*4,)
    
    '''
    
    s_t = (mtimes(x_t,warr) + mtimes(h_tm1,uarr) + barr.reshape((1, barr.shape[0])))
    hunit = uarr.shape[0]
    i  = sigmoid(s_t[:,:hunit])
    f  = sigmoid(s_t[:,1*hunit:2*hunit])
    # _c = np.tanh(s_t[:,2*hunit:3*hunit])
    _c = sigmoid(s_t[:,2*hunit:3*hunit])
    o  = sigmoid(s_t[:,3*hunit:])
    c_t = i*_c + f*c_tm1
    # h_t = o*np.tanh(c_t)
    h_t = o*sigmoid(c_t)
    return h_t, c_t

def LSTMlayer_tanh(warr, uarr, barr, x_t,h_tm1,c_tm1):
    
    '''
    c_tm1 = np.array([0,0]).reshape(1,2)
    h_tm1 = np.array([0,0]).reshape(1,2)
    x_t   = np.array([1]).reshape(1,1)
    
    warr.shape = (nfeature,hunits*4)
    uarr.shape = (hunits,hunits*4)
    barr.shape = (hunits*4,)
    
    '''
    
    s_t = (mtimes(x_t,warr) + mtimes(h_tm1,uarr) + barr.reshape((1, barr.shape[0])))
    hunit = uarr.shape[0]
    i  = sigmoid(s_t[:,:hunit])
    f  = sigmoid(s_t[:,1*hunit:2*hunit])
    _c = np.tanh(s_t[:,2*hunit:3*hunit])
    # _c = sigmoid(s_t[:,2*hunit:3*hunit])
    o  = sigmoid(s_t[:,3*hunit:])
    c_t = i*_c + f*c_tm1
    h_t = o*np.tanh(c_t)
    # h_t = o*sigmoid(c_t)
    return h_t, c_t


# %% Define Model 1
model_1 = tf.keras.models.load_model("Modeling/p20_f1_10s_sig_2h_012_027_5_150_e8_n")
def M1(x):
    weightLSTM_1 = model_1.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightD = model_1.layers[1].get_weights()
    kernel2, bias2 = weightD
    
    n1 = 10
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias2.shape[0]):
        temp2 = mtimes(h_t1[-1,:], kernel2)[i]
    temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
    out2 = sigmoid(temp2)     
    
    return out2
    
# %% Define sub model 2
model_2 = tf.keras.models.load_model('Modeling/p20_f1_10s_tanh_2h_021_032_8_90_e8_n')

def M2(x):
    weightLSTM_1 = model_2.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightD = model_2.layers[1].get_weights()
    kernel2, bias2 = weightD
    
    n1 = 10
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias2.shape[0]):
        temp2 = mtimes(h_t1[-1,:], kernel2)[i]
    temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
    out2 = tanh(temp2)     
    
    return out2

# %% Define sub model 3
model_3 = tf.keras.models.load_model("Modeling/p20_f1_5s_3s_sig_2h_029_038_4_35_e7_n")

def M3(x):
    weightLSTM_1 = model_3.layers[0].get_weights()
    warr1, uarr1, barr1 = weightLSTM_1

    weightLSTM_2 = model_3.layers[1].get_weights()
    warr2, uarr2, barr2 = weightLSTM_2
    
    weightD = model_3.layers[2].get_weights()
    kernel3, bias3 = weightD
    
    n1 = 5
    # Initialise the model with the hidden states
    h_t1 = SX.zeros((1, n1))
    # Long term memory
    C_t1 = SX.zeros((1, n1))
    for i in range(x.shape[0]):
        xi = x[i,:].reshape((1,2))
        L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
        h_t1 = vertcat(h_t1, L1[0])
        C_t1 = vertcat(C_t1, L1[1])
    
    n2 = 3
    # Initialise the model with the hidden states
    h_t2 = SX.zeros((1, n2))
    # Long term memory
    C_t2 = SX.zeros((1, n2))
    for i in range(h_t1.shape[0]):
        xi = h_t1[i,:].reshape((1,n1))
        L2 = LSTMlayer_sig(warr2, uarr2, barr2, xi, h_t2[i,:], C_t2[i,:])
        h_t2 = vertcat(h_t2, L2[0])
        C_t2 = vertcat(C_t2, L2[1])
    #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
    # Dense layer
    out2 = []
    for i in range(bias3.shape[0]):
        temp2 = mtimes(h_t2[-1,:], kernel3)[i]
    temp2 = bias3.reshape((1, bias3.shape[0])) + mtimes(h_t2[-1,:], kernel3)
    out2 = sigmoid(temp2)
    
    return out2



# # %% Define Model 1
# model_1 = tf.keras.models.load_model('p20_f1_10s_tanh_2h_013_027_3.05_1.94_e7')
# def M1(x):
#     weightLSTM_1 = model_1.layers[0].get_weights()
#     warr1, uarr1, barr1 = weightLSTM_1

#     weightD = model_1.layers[1].get_weights()
#     kernel2, bias2 = weightD
    
#     n1 = 2
#     # Initialise the model with the hidden states
#     h_t1 = SX.zeros((1, n1))
#     # Long term memory
#     C_t1 = SX.zeros((1, n1))
#     for i in range(x.shape[0]):
#         xi = x[i,:].reshape((1,2))
#         L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
#         h_t1 = vertcat(h_t1, L1[0])
#         C_t1 = vertcat(C_t1, L1[1])
    
#     #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
#     # Dense layer
#     # out2 = []
#     # for i in range(bias2.shape[0]):
#     #     temp2 = mtimes(h_t1[-1,:], kernel2)[i]
#     temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
#     out2 = tanh(temp2)     
    
#     return out2
    
# # %% Define sub model 2
# model_2 = tf.keras.models.load_model('p20_f1_10s_tanh_2h_021_030_5.56_63.9_e8')

# def M2(x):
#     weightLSTM_1 = model_2.layers[0].get_weights()
#     warr1, uarr1, barr1 = weightLSTM_1

#     weightD = model_2.layers[1].get_weights()
#     kernel2, bias2 = weightD
    
#     n1 = 10
#     # Initialise the model with the hidden states
#     h_t1 = SX.zeros((1, n1))
#     # Long term memory
#     C_t1 = SX.zeros((1, n1))
#     for i in range(x.shape[0]):
#         xi = x[i,:].reshape((1,2))
#         L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
#         h_t1 = vertcat(h_t1, L1[0])
#         C_t1 = vertcat(C_t1, L1[1])
    
#     #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
#     # Dense layer
#     # out2 = []
#     # for i in range(bias2.shape[0]):
#     #     temp2 = mtimes(h_t1[-1,:], kernel2)[i]
#     temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
#     out2 = tanh(temp2)     
    
#     return out2

# # %% Define sub model 3
# model_3 = tf.keras.models.load_model('p20_f1_10s_sig_2h_028_038_1.11_4.17_e6')

# def M3(x):
#     weightLSTM_1 = model_3.layers[0].get_weights()
#     warr1, uarr1, barr1 = weightLSTM_1

#     weightD = model_3.layers[1].get_weights()
#     kernel2, bias2 = weightD
    
#     n1 = 10
#     # Initialise the model with the hidden states
#     h_t1 = SX.zeros((1, n1))
#     # Long term memory
#     C_t1 = SX.zeros((1, n1))
#     for i in range(x.shape[0]):
#         xi = x[i,:].reshape((1,2))
#         L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
#         h_t1 = vertcat(h_t1, L1[0])
#         C_t1 = vertcat(C_t1, L1[1])
    
#     #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
#     # Dense layer
#     # out2 = []
#     # for i in range(bias2.shape[0]):
#     #     temp2 = mtimes(h_t1[-1,:], kernel2)[i]
#     temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
#     out2 = sigmoid(temp2)     
    
#     return out2
# %% Rearrange NN data

def multivariate_data(dataset, target, start_index, end_index, history_size,
                      target_size,u_shape): #, single_step=False):
    data = []
    labels = []

    start_index = start_index + history_size
    if end_index is None:
        end_index = len(dataset) - target_size
        
    for i in range(start_index, end_index):
        temp_p = dataset[i-history_size:i,:]
        for j in range(target_size-1):
        # temp_p = np.concatenate(((dataset[i-history_size:i,:]).flatten(), 
                                  # (dataset[i:i+target_size-1,:u_shape]).flatten()))
            temp = np.array([dataset[i+j,0], dataset[i,1]]).reshape((1,2))
            temp_p = np.concatenate((temp_p,temp))
        data.append(temp_p)
        labels.append((target[i:i+target_size].flatten()))
    

    return np.array( data), np.array(labels)
# %% Define dataset for 2nd layer NN
y_index = 19
past = 20
future = 1
n_model = 3
# u = np.loadtxt('u_NN_train_crop_2hr_impulse_f30_small.txt')[:120000]
# y = np.loadtxt('y_NN_train_crop_2hr_impulse_f30_small.txt')[:120001,y_index]
# uv = np.loadtxt('u_NN_train_crop_2hr_impulse_f30_small.txt')[120000:]
# yv = np.loadtxt('y_NN_train_crop_2hr_impulse_f30_small.txt')[120000:,y_index]


u = np.loadtxt("u_30hrLF_30000.txt")[:]
y = np.loadtxt("y_30hrLF_30000.txt")[:,y_index]
uv = np.loadtxt("u_vali_80hrLF_2layered_1000_010_034.txt")
yv = np.loadtxt("y_vali_80hrLF_2layered_1000_010_034.txt")[:, y_index]
# uv = np.loadtxt("u_vali_20hrLF_2layered_1000_012_037.txt")
# yv = np.loadtxt("y_vali_20hrLF_2layered_1000_012_037.txt")[:,y_index]
# Test the performance with data for M3
# u = np.loadtxt('u_NN_train_crop_2hr_impulse_f100_LPV.txt')[10100:28000]
# y = np.loadtxt("y_NN_train_crop_2hr_impulse_f100_LPV.txt")[10100:28001,y_index]
# uv = np.loadtxt('u_NN_train_crop_2hr_impulse_f100_LPV.txt')[28000:30000]
# yv = np.loadtxt("y_NN_train_crop_2hr_impulse_f100_LPV.txt")[28000:30001,y_index]



# soilPars = loamySoil()  # The soil parameters
# totalDepth, axialNodes, axialDistance, totalNodes, numOfActuators, interPoints = par_RNN.spatialVariables_1D()  #Spatial parameters
# samplingTime, samplingTimeInternal, internalTimeSteps, totTimeSteps = \
#                                  par_RNN.temporalVariable_dataGen(25000*6*60)
np.random.seed(5)

yn_train = 0.1*(max(y) - min(y))*np.random.choice([-1,1],(y.shape))*np.random.random((y.shape))
y_n = y*(1+yn_train)

y_train = y[1:]
x_train = np.zeros((u.shape[0],2))
x_train[:,0] = u
x_train[:,1] = y_n[:-1]

scaler_x = MinMaxScaler()
x_train = scaler_x.fit_transform(x_train)
scaler_y = MinMaxScaler()
y_train = scaler_y.fit_transform(y_train.reshape((y_train.shape[0],1)))

x_data, y_data = multivariate_data(x_train, y_train, 0, None, past, future,1)
train = np.zeros((x_data.shape[0], n_model))
train[:,0] = model_1.predict(x_data).reshape(x_data.shape[0],)
train[:,1] = model_2.predict(x_data).reshape(x_data.shape[0],)
train[:,2] = model_3.predict(x_data).reshape(x_data.shape[0],)
# train[:,3] = x_data[:,-1,1]


yn = 0.2*(max(yv) - min(yv))*np.random.choice([-1,1],(yv.shape))*np.random.random((yv.shape))
yv_n = yv*(1+yn)

y_test = yv[1:]
x_test = np.zeros((uv.shape[0],2))
x_test[:,0] = uv
x_test[:,1] = yv_n[:-1]


scaler_xt = MinMaxScaler()
x_test = scaler_xt.fit_transform(x_test)
y_test = y_test.reshape((y_test.shape[0],1))
scaler_yt = MinMaxScaler()
y_test = scaler_yt.fit_transform(y_test)

x_t, y_t = multivariate_data(x_test, y_test, 0, None, past, future,1)
vali =  np.zeros((x_t.shape[0], n_model))
vali[:,0] = model_1.predict(x_t).reshape(x_t.shape[0],)
vali[:,1] = model_2.predict(x_t).reshape(x_t.shape[0],)
vali[:,2] = model_3.predict(x_t).reshape(x_t.shape[0],)
vali[:,3] = x_t[:,-1,1]
# train[:,3] = x_data[:,-1,0]
# train[:,3] = model_4.predict(x_data).reshape(x_data.shape[0],)
# train[:,4] = model_5.predict(x_data).reshape(x_data.shape[0],)
# vali[:,0] = scaler_yt.transform(model_1.predict(x_t)*(0.27 - 0.13)+0.13).reshape(x_t.shape[0],)
# vali[:,1] = scaler_yt.transform(model_2.predict(x_t)*(0.32-0.21)+0.21).reshape(x_t.shape[0],)
# vali[:,2] = scaler_yt.transform(model_3.predict(x_t)*(0.38-0.28) + 0.28).reshape(x_t.shape[0],)

# vali[:,3] = x_t[:,-1,0]

# vali[:,3] = model_4.predict(x_t).reshape(x_t.shape[0],)
# vali[:,4] = model_5.predict(x_t).reshape(x_t.shape[0],)
# def alpha(x):
#     w = 
#     a = 
# class custom_loss(tf.keras.losses.Loss):
    
#     def __init__(self, y_sub = tf.zeros(5)):
#         self.y_sub = y_sub
        
#     def call(self, y_true, y_pred):
        
#         y_est = tf.linalg.matmul(y_pred, self.y_sub)
#         custom_loss=kb.square(y_actual-y_est)
#         return custom_loss

# %%
OUTPUT_SIZE = 1
model = tf.keras.Sequential() 

model.add(tf.keras.layers.Dense(3, input_shape=(3,), activation = "sigmoid"))  #, activation = 'sigmoid'

# add output layer
model.add(tf.keras.layers.Dense(3, activation = "sigmoid"))
model.add(tf.keras.layers.Dense(5, activation = "sigmoid"))
# model.add(tf.keras.layers.Dense(4, activation = "sigmoid"))
# model.add(tf.keras.layers.Dense(4, activation = "sigmoid"))
model.add(tf.keras.layers.Dense(OUTPUT_SIZE, activation = "sigmoid"))

opt = tf.keras.optimizers.Nadam(learning_rate=0.005)
model.compile(optimizer=opt, loss="mse", metrics=['mean_squared_error'])

epoch = 500

callback = tf.keras.callbacks.EarlyStopping(monitor='loss', patience=5)
Dense_history = model.fit(train,y_data.reshape((y_data.shape[0],y_data.shape[1])), epochs=epoch,
                                batch_size=200,
                                validation_split=0.3,
                            # steps_per_epoch = steps_per_epoch,
                                # validation_data = (x_t, y_t),
                                # validation_steps = validation_steps,
                                callbacks = [callback])
# for i in range(u_train.shape[0])： 
# model.save("2"2nd_layer_2hr_3s_5s_5t_1s_2"nd_layer_2hr_3s_5s_3t_2s_1s")
# %% Single Step Ahead prediction performance
# model = tf.keras.models.load_model("2nd_layer_2hr_3s_5s_5t_1s")
# model_SNN = tf.keras.models.load_model('p20_f1_10s_sig_2h_028_038_1.11_4.17_e6')
# model = tf.keras.models.load_model("Modeling/2nd_layer_4in_5s_s_10nt")
# y_pt_1 = model.predict(train[:1000,:])
y_pt_1 = model.predict(vali)
# y_pt_SNN = model.predict(x_t)
ypt_1 = scaler_yt.inverse_transform(y_pt_1)
# ypt_2 = scaler_yt.inverse_transform(y_pt_2)
# ypt = (ypt_1 + ypt_2)/2
# y_act = scaler_yt.inverse_transform(y_data[:1000,:].reshape((1000, y_t.shape[1])))   #
y_act = scaler_yt.inverse_transform(y_t.reshape((y_t.shape[0], y_t.shape[1])))
csfont = {'fontname':'Times New Roman'}
font = font_manager.FontProperties(family='Times New Roman',
                                    # weight='bold',
                                    style='normal', size=17)
plt.figure(2)
plt.plot(y_act[:2000,0], label = 'Actual Data')
plt.plot(ypt_1[:2000,0], label = 'Single-step-ahead NN')
# plt.plot(ypt_1[-500:,-1], label = 'Multi-step-ahead NN')
# plt.plot(ypt_2[30:,0], label = 'Multi-step-ahead NN -1st')
# plt.plot(ypt_2[:,-1], label = 'Multi-step-ahead NN - 30th')
plt.xlabel('Time Steps', **csfont, fontsize = 17)
plt.ylabel('Soil Moisture ($m^3/m^3$)', **csfont, fontsize = 17)
plt.legend(frameon=False, ncol=1, prop=font)

error = (sum((ypt_1 - y_act)**2/y_act**2))/ypt_1.shape[0]
print(error)    
    
# %%
# model_a = model
def Ma(x):
    weightD1 = model.layers[0].get_weights()
    kernel1, bias1 = weightD1

    weightD2 = model.layers[1].get_weights()
    kernel2, bias2 = weightD2
    
    weightD3 = model.layers[2].get_weights()
    kernel3, bias3 = weightD3
    
    weightD4 = model.layers[3].get_weights()
    kernel4, bias4 = weightD4
    
    # weightD5 = model.layers[4].get_weights()
    # kernel5, bias5 = weightD5
    
    temp1 = bias1.reshape((1, bias1.shape[0])) + mtimes(x.reshape((1,x.shape[0])), kernel1)
    out1 = sigmoid(temp1)   
   
    temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(out1, kernel2)
    out2 = sigmoid(temp2)     
    
    temp3 = bias3.reshape((1, bias3.shape[0])) + mtimes(out2, kernel3)
    out3 = sigmoid(temp3)  
    
    temp4 = bias4.reshape((1, bias4.shape[0])) + mtimes(out3, kernel4)
    out4 = sigmoid(temp4)  
    
    # temp5 = bias5.reshape((1, bias5.shape[0])) + mtimes(out4, kernel5)
    # out5 = sigmoid(temp5)  
    
    return out4
# %% Multi-step ahead prediction performance
# model_SNN = tf.keras.models.load_model('p20_f1_10s_sig_2h_028_038_1.11_4.17_e6')
# model = tf.keras.models.load_model("2nd_layer_2hr_3s_5s_5t_1s_2")
a = 900
t_pred = 20  # 30 sampling times 
error = np.zeros(a)
y_save = np.zeros((a, t_pred))

for t0 in range(43, a+43):
    print(t0)
    y_sim_1 = np.zeros((t_pred,1))
    NN0_1 = x_t[t0,:]
    y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1,1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
    y_sim_1[0,:] = np.array(Ma(y_sub)).flatten()
    u_sim = x_t[t0:t0+t_pred,-1,0]
    for i in range(t_pred-1):
        NN0_1 = np.concatenate((NN0_1[1:], np.array([u_sim[i], y_sim_1[i,0]]).reshape((1,2))))
        y_sub = np.array([DM(M1(NN0_1)), DM(M2(NN0_1)), DM(M3(NN0_1)), NN0_1[-1,1]])  #, DM(M4(NN0_1)), DM(M5(NN0_1))])
        y_sim_1[i+1,:] = np.array(Ma(y_sub)).flatten()
    yp_1 = scaler_yt.inverse_transform(y_sim_1)
    y_save[t0-43,:] = yp_1.flatten()
    error[t0-43] = sum((y_act[t0:t0+t_pred] - yp_1)**2/y_act[t0:t0+t_pred])
# yp_2 = scaler_yt.inverse_transform(y_sim_2)
# yp = (yp_1 + yp_2)/2

font = font_manager.FontProperties(family='Times New Roman',
                                   # weight='bold',
                                    style='normal', size=15)
a = 900
tplot = np.arange(a)
plt.figure(1)
plt.plot(tplot[:], y_act[:a], label = "Actual Trajectory", linewidth = 4, color = "k")
plt.plot(tplot[:], y_save[:a,0], label = "Single-step-ahead Prediction")
plt.plot(tplot[9:], y_save[:a-9,4], label = "Five-step-ahead Prediction", color = "y")
plt.plot(tplot[19:], y_save[:a-19,19], label = "Twenty-step-ahead Prediction", color = "g")
plt.legend(frameon=False, ncol=1, prop=font)
plt.legend(frameon=False, ncol=1, prop=font)
plt.xlabel("Time Steps", **csfont, fontsize = 17)
plt.ylabel('Soil Moisture ($m^3/m^3$)', **csfont, fontsize = 17)

# e_0 = np.sqrt(sum((y_act[:900,0] - y_save[:900,0])**2/y_act[:900,0]**2))/900 # 0.0027434248213169844
# e_5 = np.sqrt(sum((y_act[4:904,0] - y_save[:900,4])**2/y_act[4:904,0]**2))/900  # 0.005014430666105703
# e_10 = np.sqrt(sum((y_act[9:909,0] - y_save[:900,9])**2/y_act[9:909,0]**2))/900 # 0.00700915468836735
# e_20 = np.sqrt(sum((y_act[19:919,0] - y_save[:900,19])**2/y_act[19:919,0]**2))/900 # 0.008967830204436932

# np.savetxt("Modeling/y_pred_80hrLF_2layered_1000_010_034_20nv_2layer_sub_noise.txt",y_save)
# model.save("Modeling/2nd_layer_2hr_3s_5s_5t_1s_noise")
e_0 = np.sqrt(sum((y_act[:900,0] - y_save[:900,0])**2)/900)/(max(y_act[:900,0] - y_save[:900,0]))
e_5 = np.sqrt(sum((y_act[4:904,0] - y_save[:900,4])**2)/900)/(max(y_act[:900,0] - y_save[:900,0])) 
e_10 = np.sqrt(sum((y_act[9:909,0] - y_save[:900,9])**2)/900)/(max(y_act[:900,0] - y_save[:900,0]))
e_20 = np.sqrt(sum((y_act[19:919,0] - y_save[:900,19])**2)/900)/(max(y_act[:900,0] - y_save[:900,0]))
# t0=552
# plt.figure(5)
# plt.plot(y_act[t0:t0+t_pred,-1], label = 'Actual Data', color = "k")  # , marker = 'o'
# plt.plot(ypt_1[t0:t0+t_pred,-1], label = 'Single-step', color = "b")  # , marker = 'd'
# plt.plot(yp_1[0:,-1], label = 'Multi-step', color = "r")  # , marker = '^'
# plt.xlabel('Time Steps', **csfont, fontsize = 17)
# plt.ylabel('Soil Moisture ($m^3/m^3$)', **csfont, fontsize = 17)

# plt.legend(frameon=False, ncol=1, prop=font)

# %%
# class explicitNN:
#     def __init__(self, model):
#         self.x = np.zeros((20,2))
#         self.model = model
#     def explicit(self):
#         weightLSTM_1 = self.model.layers[0].get_weights()
#         warr1, uarr1, barr1 = weightLSTM_1

#         weightD = self.model.layers[1].get_weights()
#         kernel2, bias2 = weightD
    
#         n1 = 10
#         # Initialise the model with the hidden states
#         h_t1 = SX.zeros((1, n1))
#         # Long term memory
#         C_t1 = SX.zeros((1, n1))
#         for i in range(self.x.shape[0]):
#             xi = self.x[i,:].reshape((1,2))
#             L1 = LSTMlayer_sig(warr1, uarr1, barr1, xi, h_t1[i,:], C_t1[i,:])
#             h_t1 = vertcat(h_t1, L1[0])
#             C_t1 = vertcat(C_t1, L1[1])
    
#     #  array([-2.7764497, -2.6578658, -2.6829786, -2.7356846, -2.73696  ]
#         # Dense layer
#         out2 = []
#         for i in range(bias2.shape[0]):
#             temp2 = mtimes(h_t1[-1,:], kernel2)[i]
#             temp2 = bias2.reshape((1, bias2.shape[0])) + mtimes(h_t1[-1,:], kernel2)
#             out2 = tanh(temp2)     
#         return out2  
    

# M1 = explicitNN(model_1)
# M2 = explicitNN(model_2)
# M3 = explicitNN(model_3)
# M4 = explicitNN(model_4)
# M5 = explicitNN(model_5)